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Introduction 

It  is  well  documented  from  laboratory  experiments  that  the  attenuation  of  elastic  waves 
is  amplitude  dependent  at  strains  greater  than  I  O'6.  There  exists  therefore  an  intermediate 
strain  regime  in  which  the  deformation  is  not  strong  enough  to  cause  macroscopic  material 
failure,  but  is  large  enough  to  affect  wave  propagation.  Because  the  rheology  depends  on 
amplitude,  wave  propagation  is  nonlinear,  and  reflects  the  influence  of  inelastic 
mechanisms.  In  the  high-strain  regime,  these  mechanisms  dominate  and  control  the  rate  of 
energy  dissipation,  but  in  the  low-strain  regime,  their  effects  are  masked  by  linear, 
amplitude  independent  anelastic  mechanisms. 

The  nonlinearity  means  that  the  analytical  tools  usually  employed  to  solve  wave 
propagation  problems  are  not  valid  in  this  intermediate-strain  regime.  This  raises  a  number 
of  questions  regarding  the  acquisition  and  interpretation  of  experimental  data.  For  instance 

•  In  resonating  bar  experiments,  Lorentz  peaks  become  quite  distorted  and 
asymmetric,  and  the  interpretation  of  the  half-power  peak  width  in  terms  of  Q'1  is 
not  strictly  valid,  since  the  predicted  width  of  the  Lorentz  peak  stems  from  a  linear 
perturbation  treatment  of  the  complex  Helmholtz  equation. 

•  In  pulse  propagation  experiments,  the  use  of  spectral  ratios  to  estimate  Q'1  depends 
implicitly  on  the  validity  of  a  linear  theory.  The  use  of  rise  times  as  Q'1  estimators 
is  also  open  to  question  since  the  relationship  between  rise  time  and  attenuation  has 
been  derived  only  in  the  linear  case. 

Nevertheless,  linear  analysis  of  laboratory  data  which  clearly  exhibit  nonlinear  traits  is 
attractive,  in  the  sense  that  it  yields  fairly  simple  and  apparently  self-consistent  results. 
Furthermore,  a  reasonable  physical  interpretation  appears  feasible:  their  own  experiments 
have  led  Stewart  et  al.  [1983]  to  propose  a  physical  model  based  on  Hertzian  contacts 
across  cracks  which  generalizes  the  earlier  models  of  Walsh  [1966]  and  Mavko  [1979]. 

The  nonlinear  regime  must  prevail  within  a  range  of  scaled  distances  from  a  contained 
explosive  source,  since  strain  amplitudes  in  the  very  near  field  are  sufficient  to  cause 
material  failure,  while  far-field  radiation  falls  clearly  in  the  linear  low-amplitude  regime. 
This  intermediate  strain,  nonlinear  zone  contains  the  transition  from  what  has  traditionally 
been  called  the  “elastic”  radius  in  numerical  source  models  to  the  far-field,  low-strain 
(an)elastic  regime.  Minster  and  Day  [1986]  attempted  to  examine  this  and  other  issues  by 
comparing  the  decay  of  peak  displacements  and  velocities  for  the  Cowboy  series  with 
laboratory  results  on  salt  obtained  by  Tittman  [1983].  However,  the  algorithm  they  used  is 
inadequate  since  it  does  not  satisfy  local  conditions  both  in  time  and  space,  and  thus  is 
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nonphysical.  A  better  algorithm  can  be  devised  using  the  Pad6  approximations  developed 
by  Day  and  Minster  [1984, 1990],  or  the  least-squares  approximations  of  Witte  [1989]  and 
Witte  and  Richards,  [1990],  although  any  procedure  which  depends  on  selecting  Q  as  a 
particular  measure  of  (amplitude-dependent)  attenuation  is  likely  to  be  an  inadequate 
representation  of  the  true  physical  processes. 

In  this  paper,  we  present  a  series  of  simple,  one-dimensional  numerical  experiments 
designed  to  test  the  consistency  of  the  models  and  laboratory  results.  When  it  comes  to 
wave  propagation  in  a  nonlinear  medium,  intuition  developed  from  linear  Q  models  is  a 
rather  poor,  often  misleading  guide,  and  our  results  are  mixed,  and  do  not  seem  to  be  easily 
predictable  from  simple  arguments.  For  example,  Q'1  estimates  derived  from  the  half¬ 
width  of  a  resonance  peak  appear  to  be  surprisingly  accurate  well  into  the  nonlinear  regime. 
Similarly,  the  nonlinear  interaction  of  one-dimensional  pulses  does  not  lead  to  strong 
departures  from  linear  superposition  in  the  range  of  nonlinear  behavior  we  have 
considered.  On  the  other  hand,  the  propagation  of  a  one-dimensional  narrow  pulse  through 
a  medium  with  frequency-independent,  but  amplitude-dependent,  Q  is  not  described 
accurately  by  an  equivalent  Q  operator. 

Laboratory  constraints 

Seismic  attenuation  depends  on  a  large  number  of  conditions.  These  include  pressure, 
temperature,  porosity  and  saturation  of  the  rock,  and  of  course  frequency  and  amplitude  of 
the  propagating  wave.  The  bulk  of  laboratory  evidence  accumulated  to  date  indicates  that 
nonlinear  processes  are  activated  in  rocks  at  strain  levels  above  10'6  [e.g.  Mavko,  1979; 
Johnston  and  Toksdz,  1980;  Liu  and  Peselnick,  1979;  Bulau  et  al.,  1984],  Stewart  et  al. 
[1983]  reviewed  the  available  data  and  conducted  additional  experiments  in  which 
nonlinearity  was  shown  to  decrease  markedly  at  confining  pressures  high  enough  to  close 
many  cracks.  They  explain  the  observations  with  a  model  where  attenuation  for  dry  rock  at 
low  to  moderate  confining  pressure  is  caused  by  frictional  work  dissipation  at  crack 
asperities  in  Hertzian  contact.  For  large  strains,  the  attenuation  increases  linearly  with 
strain  amplitude  and  crack  density,  and  decreases  with  confining  pressure  according  to  the 
relation: 


Q'l=Qa  +k£P4fie 


(1) 


where  k  is  a  material  constant,  £  is  the  crack  density  ,  e  is  the  strain  amplitude,  and  P  the 
confining  pressure.  Here  Qa  represents  the  linear  attenuation  controlled  by  anelastic 
mechanisms  that  mask  the  nonlinear  ones  at  low  strain.  Tittman  [1983]  showed  that  a 
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similar  amplitude  dependence  also  holds  for  halite,  both  in  the  form  of  pressed  salt  and  in 
dome  salt. 

The  experimental  and  analytical  methods  used  to  measure  Q  in  the  laboratory  fall 
generally  in  three  main  categories  fe.g.  Nowick  and  Berry,  1972J: 

1 .  Quasi-static  measurements,  in  which  we  can  include  hysteresis  loop  techniques, 
stress  relaxation  and  creep  recovery  experiments. 

2.  Resonance  experiments,  in  which  a  device  with  large  external  inertia  is  used  to 
control  the  period  of  free  or  forced  resonance  of  the  system.  This  class  includes 
torsion  pendulum  measurements  and  vibrating  bar  measurements. 

3.  Wave  propagation  experiments,  which  include  both  continuous  wave  and  short 
pulse  measurements. 

A  important  issue  for  our  present  purposes  is  that  the  experimental  setup  and  data 
reduction  techniques  used  in  many  of  these  experiments  are  based  on  linear  concepts.  For 
instance,  in  one  type  of  hysteresis  loop  measurements,  the  phase  angle  <p  of  stress-strain 
cycles  is  evaluated  by  signal  analysis,  and  attenuation  is  measured  by  Q~l  =  tan  (f>  which 
holds  as  long  as  the  rheology  is  linear.  Similarly,  in  forced  resonance  experiments,  the  half 
width  Aft)  of  the  “Lorentz”  peak  is  used  to  define  Q1  at  or  near  the  resonance  frequency 
(Of),  by  the  relation  Q'^Oo)  =  Aft)  /  (Do-  Again,  this  equation  is  derived  by  linear 
perturbation  analysis  of  the  resonant  system,  and  holds  in  the  low-loss  case.  In  wave 
propagation  experiments,  typical  approaches  involve  the  calculation  of  spectral  ratios  of 
pulses  propagated  through  the  medium  under  study  and  through  a  high-Q  reference 
medium.  Another  approach  is  to  measure  the  increase  in  rise-time  of  a  pulse  as  it  is 
attenuated.  Both  methods  rely  implicitly  on  a  linear  analysis  of  pulse  propagation.  Of 
concern  is  the  fact  that  most  documented  observations  of  amplitude-dependent  Q  rely  on 
such  techniques,  and  the  relation  (1),  which  clearly  indicates  nonlinear  behavior,  is  actually 
fitted  to  observations  derived  from  a  linear  analysis.  It  is  therefore  legitimate  to  raise  the 
question  of  whether  the  interpretation  of  amplitude-dependent  Q  observations  contains 
internal  inconsistencies,  and  possible  biases.  Furthermore,  whereas  our  intuition  is  usually 
a  very  good  guide  when  it  comes  to  gauging  the  effects  of  small  perturbations  in  the  linear 
regime,  this  is  no  longer  the  case  when  nonlinearities  dominate.  Consequently,  it  is  not  a 
foregone  conclusion  that  the  effects  of  nonlinear  Q  can  be  gauged  with  “back-of-the- 
envelope”  arguments. 

A  notable  exception  to  the  concerns  raised  above  consists  of  hysteresis  loop 
experiments  in  which  attenuation  is  defined  by  Q'1  =  AE  /  4n<E>  [e.g.  O'Connell  and 
Budiansky  ,  1978],  where  AE  is  the  dissipated  energy  (measured  by  the  loop  area)  and 
<E>  the  average  stored  energy  in  the  course  of  a  cycle.  The  definition  of  Q’1  does  not 
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depend  on  the  shape  of  the  loop,  which  is  elliptical  in  the  linear  regime  and  develops  cusps 
when  nonlinearity  appears  at  large  strains  [e.g.  Brennan  and  Stacey,  1977;  Liu  and 
Peselnick ,  1983].  Coyner  and  Martin  [1990]  have  recently  shown  that  a  correct 
interpretation  of  such  hysteresis  loops  leads  to  Q~]  estimates  for  a  variety  of  rock  types 
which  are  quite  consistent  with  estimates  obtained  by  other  techniques.  This  includes  the 
monotonic  increase  of  Q'1  with  strain  amplitude  at  strains  greater  than  10~6.  This 
consistency  may  be  taken  as  a  justification  to  accept  published  studies  on  nonlinear  Q  in 
spite  of  the  use  of  linear  data  analysis. 

In  an  attempt  to  identify  possible  pitfalls  and  to  evaluate  quantitatively  the  predicted 
effects  of  nonlinearity,  we  have  conducted  three  groups  of  simple  numerical  experiments 
with  one-dimensional  linear  and  nonlinear  systems  of  the  type  that  can  ber  described  by 
equation  (1).  They  are: 

1 .  simulations  of  forced  vibration  of  nonlinearly  damped  one-dimensional  oscillators, 
to  examine  the  distortion  of  Lorentzian  peaks  induced  by  nonlinearity, 

2.  tests  of  whether  a  pair  of  one-dimensional  pulses  propagating  in  opposite  directions 
in  a  nonlinear  medium  superpose  quasi-linearly, 

3.  comparison  of  the  results  of  spectral  ratio  analysis  and  pulse  shape  interpretations 
(e.g.  rise  times)  in  both  linear  and  nonlinear  one-dimensional  pulse  propagation 
simulations. 

Numerical  experiments 

Model  formulation 

To  perform  the  suite  of  experiments  described  above,  we  use  an  attenuation  model 
described  by 

Q-l=Q-l+ye  (  (2) 

where  e  is  the  strain  amplitude  and  yis  a  material  constant.  This  is  a  specialization  of 
equation  (1),  which  was  introduced  by  Mavko  [1979]  and  was  subsequently  used  by 
Minster  and  Day  [  1986]. 

For  simulations  of  resonance  experiments,  we  adopt  a  3-parameter  attenuation  model 
characterized  by  a  single  relaxation  mechanism  for  simplicity.  This  mechanism  is  selected 
to  center  the  absorption  band  at  the  resonance  frequency,  so  that  the  appropriate  Q'1  value 
is  that  at  the  center  of  the  absorption  band.  In  that  case  the  attenuation  model  is  completely 
specified  by  coo,  the  center  frequency  of  the  band,  Qo ~  Q~!(coo),  the  attenuation  at  the 
center  of  the  band,  and  by  a  modulus  or  a  wave  velocity.  In  terms  of  a  standard  linear  solid 
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model  [e.g.  Nowick  and  Berry,  1972;  Liu  et  al.,  1976],  the  strain  retardation  and  stress 
relaxation  times  are  given  by: 


Tf  =  o^[^Qq+  1  +  Qo) 

(3) 

ra=  %{Vq02+  1  -  Qo) 

(4) 

Since  we  are  interested  in  solving  a  quasi  monochromatic  problem,  we  can  specify  the 
phase  velocity  cq  at  the  center  frequency  cup  as  the  third  model  parameter.  The  relaxed  and 
unrelaxed  moduli  are  then  given  by: 


Mr 


2pd  .  M  = 

1  +  '  U  1  +  TcJXe 


(5) 


Although  this  model  entails  a  frequency  dependence  of  Q1,  with  Q1  increasing  as  (0  at 
frequencies  less  than  u>o ,  and  decreasing  as  (o1  at  frequencies  higher  than  coo,  the 
influence  of  frequency  dependence  is  in  fact  negligible,  because  the  width  of  the  Debye 
peak  (the  absorption  band)  is  much  larger  than  that  of  the  Lorentz  peak  (the  resonance)  in 
all  calculations  described  below.  In  other  words,  Q'1  is  essentially  frequency-independent 
in  the  range  covered  by  the  calculations.  To  simulate  nonlinearity,  we  simply  make  Qo'1 
depend  lincai ly  on  die  strain  amplitude  of  the  oscillations. 

For  pulse  propagation  experiments,  we  use  the  Pade  approximant  method  described  by 
Day  and  Minster  [1984],  slightly  modified  to  include  poles  at  zero  and  infinite  frequency, 
as  described  by  Day  and  Minster  [1990].  This  allows  us  to  use  a  broad  absorption  band 
(that  is,  broader  than  the  frequency  band  covered  by  the  finite-difference  algorithm),  and  to 
specify  the  attenuation  properties  in  each  cell  at  each  time  step  independently,  based  on  the 
current  strain  amplitude.  This  remedies  the  noncausal  character  of  the  calculations 
performed  by  Minster  and  Day  [1986]. 


Resonance  experiment 

In  resonant  bar  experiments,  amplitude  is  obtained  as  a  function  of  frequency  near  a 
normal  mode  of  the  sample.  (The  frequency  of  this  mode  is  typically  adjusted  by  addition 
of  large  external  inertia.)  The  amplitude  and  width  of  the  resulting  “Lorentzian”  peak  are 
controlled  by  the  attenuation  -  larger  attenuation  leads  to  a  lower,  broader  peak.  For  linear 
attenuation  in  which  the  width  of  the  resonance  peak  is  relatively  narrow  (e.g.  Q~l  <  ~0.1), 
the  resonance  peak  is  essentially  symmetric  (Figure  la)  and  Q  can  be  obtained  simply  from 
the  width  of  the  peak  at  the  half-power  level,  Q'1  -  (a>2-ut/)iajr.  In  the  case  of  amplitude- 
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dependent  attenuation,  the  peak  becomes  distorted  and  is  steeper  on  the  lower  frequency 
side  (see  Figure  lb).  A  simple  argument  to  explain  this  distortion  invokes  the  effect  of  the 
amplitude  dependence  on  the  phase  velocity:  the  larger  amplitude,  and  thus  larger 
attenuation,  near  the  top  of  the  peak  causes  a  reduction  in  the  phase  velocity  and  a 
corresponding  shift  of  the  resonance  to  lower  frequencies.  However  the  sides  of  the  peak 
are  not  shifted,  because  the  oscillations  have  low  amplitudes  at  those  frequencies,  so  the 
overall  peak  shape  is  asymmetric  and  skewed  to  the  left. 

Figure  2  shows  examples  of  distorted  Lorentz  peaks  obtained  from  high-strain  flexure 
experiments  in  salt  (data  from  Bulau  et  al .,  1985).  Despite  the  often  considerable 
asymmetry  in  the  peak  shapes,  experimenters  typically  continue  to  estimate  Q1  by 
measuring  the  peak  width  at  the  half-power  level.  In  order  to  test  for  possible  biases 
introduced  by  this  practice,  we  computed  synthetic  Lorentzian  peaks  through  numerical 
simulations  of  a  harmonic  oscillator  with  nonlinear  attenuation.  We  use  a  simple  finite 
difference  scheme  and,  as  described  earlier,  take  attenuation  to  be  of  the  form  (2).  When 
these  simulations  are  performed  for  constant  y,  but  increasing  forcing  amplitude,  the 
Lorentzian  peaks  become  distorted  and  are  indeed  skewed  to  the  left  (Figure  3).  However, 
when  we  calculate  Q'1  from  these  peaks  using  the  expression  Q1  =  (a)2-(0])/o}r  and 
ignoring  the  peak  asymmetry,  the  resulting  values  of  Q'1  vary  linearly  with  the  maximum 
strain  amplitude  emax  (Figure  4).  The  slope  /  computed  from  these  Q  estimates  is 
reasonably  close  to  the  /actually  used  in  the  simulations.  In  other  words,  application  of 
the  linear  concept  fields  the  "correct”  answer  to  an  acceptable  approximation.  We  conclude 
therefore  that  Q  measurements  based  on  the  half-w  idth  of  distorted  Lorentzian  peaks 
probably  give  reasonably  correct  results,  even  well  into  the  nonlinear  regime. 

In  principle,  it  might  be  possible  to  recover  the  amplitude  dependence  of  Q~l  (e.g.  )) 
from  a  single  Lorentz  peak  measurement  by  measuring  the  distortion  of  the  peak.  In  order 
to  test  this  possibility,  we  considered  a  skewness  parameier  defined  as  5  =  (0>2-0ir)Hu)r- 
a>i).  In  our  synthetic  experiments,  s  increased  with  increasing  strain  but  seemed  to  saturate 
at  about  j=/_5.  In  contrast,  the  Lorentz  peaks  from  the  laboratory  salt  experiments  (Figure 
2)  were  more  severely  distorted  with  values  of  s  up  to  2  or  higher.  Thus,  our  simulations 
were  not  completely  successful  in  reproducing  the  shapes  of  some  of  the  observed  Lorentz 
curves.  This  may  reflect  a  limitation  in  our  modeling  procedure,  or,  alternatively,  may 
indicate  that  the  attenuation  in  the  laboratory  experiments  cannot  be  modeled  as  a  simple 
linear  function  of  strain  amplitude.  In  spite  of  this  caveat,  we  conclude  that  resonant-bar 
laboratory  measurements  of  amplitude  dependence  of  Q~l  are  in  fact  a  fairly  accurate 
characterization  of  the  nonlinear  material  properties. 

(■ 


Pulse  superposition  experiment 

The  pulse  superposition  numerical  experiments  were  motivated  by  results  from  some 
small-scale  laboratory  experiments  performed  by  Larson  [1982J.  The  laboratory 
experiments  were  intended  to  test  whether  linear  superposition  applied  to  the  propagation  of 
pulses  with  maximum  strain  of  about  7  x  l(H.  At  this  strain  level,  substantial  departure 
from  linear  superposition  would  be  expected,  In  the  laboratory  experiments,  two  nearly 
identical  explosive  sources  were  simultaneously  detonated  approximately  10  cm  apart  in 
pressed  salt,  and  the  velocity  waveforms  were  recorded  along  the  plane  of  symmetry 
between  the  sources.  These  velocity  recordings  were  obtained  at  a  source-receiver  distance 
comparable  to  the  separation  distance  between  the  two  sources.  The  resulting  velocity 
records  from  the  combined  explosions  were  in  each  case  in  very  close  agreement  with  the 
sum  of  the  velocity  records  from  the  individual  charges. 

The  above  result  was  interpreted  as  evidence  for  the  applicability  of  linear  superposition 
in  wave  propagation  in  salt  at  relatively  high  strains.  This  is  in  apparent  contradiction  to 
other  laboratory  measurements  in  the  same  material,  which  consistently  show  evidence  of 
amplitude-dependent  losses  at  comparable  strain  levels.  An  alternative  interpretation, 
which  we  examine  in  this  section,  is  that  an  experiment  of  this  type  is  simply  not  sensitive 
to  deviations  from  linearity. 

To  address  this  issue,  we  describe  a  series  of  numerical  simulations  of  pulse 
propagation.  The  numerical  experiments  are  for  plane  waves  in  a  simple  nonlinear 
rheology  which  nevertheless  mimics  the  laboratory-observed  amplitude  dependence  of  Q  as 
summarized  earlier. 

The  model  for  nonlinear  attenuation  is  constructed  in  two  steps.  First,  we  use  the  Pad£ 
approximant  method  of  Day  and  Minster  [  1984,  1990)  to  convert  the  stress-strain  relation 
of  a  linear,  anelastic  solid,  with  frequency-independent  Q,  into  differential  form.  Minster 
1 1 978a |  shows  that  an  absorption  band,  with  Q  nearly  constant  at  (?o,  and  with  minimum 
and  maximum  relaxation  times  T\  and  T2,  respectively,  yields  the  following  relation 
between  stress  history,  o(t)  and  strain  history,  e(t) 
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where  Mu  is  the  unrelaxed  modulus.  Day  and  Minster  show  that  (6)  can  be  approximated 
by 


o{t)  =  Mu 


£<t)  -  £  £(/) 

i  l 


(7) 


where  the  C;  's  are  relaxation  tenns  governed  by  the  n  linear  equations 


~  *■  v,C  =  Qo  dt)  ,  (8) 

The  constants  v,  and  vv,  which  depend  on  the  order  of  approximation,  n,  are  given  by  Day 
and  Minster  [1984,  1990J,  who  also  show  that  the  operator  defined  by  (7)  and  (8) 
converges  to  the  exact  result  (6)  as  n  increases. 

The  second  step  is  to  generalize  (8)  by  introducing  a  linear  dependence  of  Qo  on  strain 
amplitude  according  to  (2): 

^  =  T|-  - ■*  vi  |gu4)i£(r)|)  £(/)  .  (9) 

Then,  (7)  and  (9)  constitute  the  stress-strain  equations  for  our  set  of  one-dimensional  finite 
difference  simulations. 

The  source  is  a  velocity  pulse  — in  agreement  with  Larson’s  experiments —  with  time 
history 


0.5  [1  -  cos{2k  tfT))] .  0<t<T  (10) 

0,  t<0  or  t>T 

introduced  at  a  plane  boundary.  Figure  5  shows  the  velocity  time  histories  at  5  different 
distances  from  the  source,  for  a  calculation  in  which  Qa  is  1000  and  the  product  of  >  and 
the  maximum  strain  at  the  source  is  equal  to  0.3.  In  other  words,  the  amplitude 
dependence  is  such  that  Q  assumes  a  minimum  value  of  about  3  near  the  source.  As  a 
result  of  the  very  low  effective  values  of  Q,  the  plane  wave  amplitude  attenuates  rapidly 
with  distance  of  propagation.  Furthermore,  the  amplitude  dependence  of  the  attenuation 
causes  the  pulse  to  evolve  into  a  distinctive  shape,  with  an  initial  steep  rise,  followed  by  a 
more  gradual  slope.  The  break  in  slope  is  not  present  for  the  linear  attenuation  model  (i.e., 
?  =  0).  As  in  the  case  of  the  distorted  Lorentz  peak.>,  we  can  offer  a  simple  explanation  for 
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this  peculiar  pulse  shape:  the  crest  of  the  pulse  has  large  amplitude  and  is  thus  associated 
with  low  Q,  and  therefore  low  phase  velocity.  Consequently,  the  crest  of  the  wave  gets 
progressively  delayed  with  respect  to  both  the  low  amplitude  onset  and  tail  of  the  pulse. 

Figure  6  shows  the  corresponding  time  histories  when  two  plane  waves  traveling  in 
opposite  directions  interact.  In  this  case,  identical  velocity  pulses  of  form  (10)  were 
introduced  simultaneously,  separated  by  20  times  the  initial  pulse  width  (travel  time 
separation  of  207).  Figure  7a  compares  the  two-source  waveform,  at  the  midpoint,  with 
twice  the  single-source  waveform.  This  is  the  one-dimensional  analogue  of  Larson's 
laboratory  superposition  experiment.  As  in  the  laboratory  experiment,  the  numerical 
experiment  shows  near-perfect  agreement.  Finally,  Figure  7b  compares  the  two-source 
waveforms,  after  pulse  interaction  is  completed,  with  the  single-source  waveform.  Again, 
agreement  is  very  close:  although  the  effect  of  imperfect  superposition  is  visible  in  the 
theoretical  calculation,  it  is  small  enough  that  it  would  be  masked  by  experimental 
uncertainties.  Comparable  results  were  obtained  in  this  comparison  when  destructive, 
rather  than  constructive,  interference  between  the  sources  was  present  in  the  numerical 
experiment. 

These  results  show  that  this  type  of  experimental  observation  is  not  diagnostic  of 
deviations  from  linearity.  This  insensitivity  apparently  results  because  the  pulses 
superpose  only  briefly  as  they  cross.  For  the  amplitude-dependent  relaxation  model,  this 
brief  interaction  is  insufficient  to  discemibly  perturb  the  individual  source  pulses.  This  is 
so  even  though  the  individual  pulses  themselves  have  been  substantially  modified  as  a 
result  of  the  nonlinearity  of  the  model. 

Spectral  ratios 

One  of  the  standard  methods  for  estimating  Q  from  seismic  observations  is  the  spectral 
ratio  method.  For  amplitude-independent  attenuation,  the  spectral  attenuation  of  a  plane 
wave  is  exponential  in  the  frequency,/,  and  propagation  distance*, 

v(fjc)  =  e y  (II) 

where  c  is  the  wave  speed.  Then,  if  we  assume  that  Q  is  approximately  frequency- 
independent,  the  logarithm  of  the  spectral  ratio  of  observations  at  distances  x\  and  *2  is 
linear  in  the  frequency /: 


(12) 
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The  Pad6  approximant  method  reproduces  this  behavior  of  frequency-and  amplitude- 
independent  Q  models,  as  shown  by  Day  and  Minster  [1990],  In  this  section,  we  examine 
the  effect  of  amplitude  dependent  Q  on  spectral  ratio  measurements  of  plane  waves. 
Numerical  experiments  are  performed  using  a  displacement  — rather  than  velocity —  source 
pulse  of  form  (10).  One  complication  which  we  examine  in  particular  results  directly  from 
the  introduction  of  nonlinearity,  and  concerns  the  possibility  of  energy  transfer  from  one 
frequency  component  to  another. 

For  comparison  purposes,  and  as  verification  of  the  numerical  method,  we  present 
results  for  a  linear  (amplitude-independent  Q )  calculation,  with  Q  equal  to  20,  in  Figures  8 
and  9.  Figure  8  shows  displacement  time  histories  at  distances  of  4,  12,  and  20  times  the 
source  pulse  width  (i.e.,  4 cT,  12 cT,  and  20 cT).  Figure  9  shows  the  corresponding  spectra 
divided  by  the  source  spectrum.  Note  that,  over  most  of  the  frequency  range  plotted,  the 
spectral  ratios  follow  the  log-linear  relationship  of  Equation  12,  as  expected.  Also  note  that 
all  spectral  ratios  approach  unity  in  the  low-frequency  limit,  to  very  high  accuracy.  This 
agrees  with  the  usual  “linear”  reasoning  that  when  the  wavelength  exceeds  the  propagation 
distance,  then  attenuation  has  little  chance  to  affect  spectral  amplitudes. 

Figure  10  shows  time -domain  displacements  for  a  nonlinear  calculation,  with  Qa  equal 
to  1000  and  the  product  [)  x  emax]  equal  to  0.3,  i.  e.,  Qmin  ~  3.  The  most  notable 
difference  between  these  time-histories  and  those  for  the  linear  case  is  that  the  tail  of  the 
pulse  is  of  shorter  duration  in  the  nonlinear  case,  making  the  pulse  more  symmetric  in 
appearance. 

The  spectral  ratios  given  in  Figure  1 1  show  a  substantial  deviation  from  the  linear  case, 
however.  The  most  important  difference  is  that  the  spectral  ratios  in  this  case  do  not 
approach  1  at  zero  frequency.  At  a  distance  of  4cT,  the  zero-frequency  spectral  value  is 
only  half  of  its  value  for  the  source  pulse,  and  the  zero-frequency  spectrum  continues  to 
diminish  by  a  further  15%  as  distance  increases  to  20c7\  Furthermore,  at  least  in  some 
cases,  the  attenuated  spectra  have  a  slight  peak,  which  means  that  the  pulse  does  not  remain 
unipolar  as  it  propagates,  a  phenomenon  which  can  be  detected  in  the  pulse  shapes  shown 
on  Figure  8.  This  low-frequency  behavior  is  a  nonlinear  phenomenon  which  has  no 
counterpart  in  the  linear  case. 

The  other  distinction  between  spectral  ratios  for  the  linear  and  nonlinear  cases  is  the 
increasing  slope  in  the  nonlinear  case  above  a  frequency  of  about  l  IT.  This  is  a 
consequence  of  the  amplitude  dependence  of  the  attenuation  mechanism.  High  frequencies 
associated  with  the  discontinuous  acceleration  at  the  onset  of  the  displacement  pulse  are 
relatively  less  attenuated  in  the  nonlinear  case,  because  the  pulse  onset  is  associated  with 
relatively  small  strain  amplitudes.  The  spectral  ratios  for  the  nonlinear  case  are  further 


10 


distorted  as  the  frequency  approaches  2/7,  as  a  result  of  energy  transfer  from  one 
frequency  band  to  another.  The  source  spectrum  has  a  null  at  2/7,  and  nonlinearity  causes 
a  small  amount  of  energy  leakage  into  this  band.  While  the  effect  is  relatively  small,  the 
presence  of  a  null  in  the  source  spectrum  results  in  a  large  effect  on  the  spectral  ratios  in  the 
vicinity  of  2/7.  As  with  the  low-frequency  behavior,  these  high-frequency  phenomena 
have  no  counterpart  in  the  case  of  linear  attenuation  models. 

For  nonlinear  simulations  with  somewhat  smaller  values  of  >,  similar  phenomena 
occur,  but  are  correspondingly  less  pronounced.  For  [?  x  emax)  equal  to  0.1  (Qmin  ~  10), 
for  example,  the  zero-frequency  spectral  ratio  is  about  0.8  at  a  distance  of  20c7. 

The  numerical  experiments  also  provide  evidence  that  the  behavior  of  the  low- 
frequency  spectral  components  is  likely  to  be  sensitive  to  the  shape  of  the  source  pulse. 
For  example,  simulations  verify  that  a  unipolar  velocity  pulse  propagates  in  the  present 
model  without  reduction  of  the  zero-frequency  spectral  component,  in  contrast  to  the  above 
results  for  a  unipolar  displacement  pulse.  In  other  words,  in  the  nonlinear  model, 
integration  does  not  commute  with  the  attenuation  operator,  as  must  be  the  case  for  all 
linearly  anelastic  models.  This  means  in  particular  that  nonlinear  pulse  propagation 
characteristics  depend  on  the  shape  (i.e.  the  spectral  content)  of  the  pulse  itself. 

Yet  another  method  for  estimating  Q1  from  pulse  propagation  experiments  is  to  map 
directly  the  rise  time  of  the  pulse  into  attenuation  estimates  through: 

Tr=  T?  +  CtuQl  (13) 

where  r,  is  the  rise  time  of  the  attenuated  pulse,  T?  the  rise  time  of  the  input  pulse,  and  tu  is 
the  (unrelaxed)  travel  time.  Here  C  is  a  constant  equal  to  approximately  0.5,  as  shown  by 
several  authors  [e.g.  Gladwin  and  Stacey,  1974;  Minster,  1978b;  Kjartansson,  1979]. 
Stewart  et  al.  (1983]  used  the  technique  in  high-strain  experiments  to  demonstrate 
amplitude  dependence  of  Q'1  in  dry  Berea  sandstone.  The  pulses  shown  in  Figures  5  and  6 
illustrate  the  problem  we  have  in  verifying  this  relation  for  the  nonlinear  case.  Because  of 
the  peculiar  distortion  of  the  pulse,  the  rise  time  is  not  easily  defined  using  the  tangent  at  the 
inflection  point  as  defined  by  Minster  [1978b].  Further,  for  linear  mechanisms,  not  only 
the  rise  time,  but  all  time  scales  within  the  pulse  scale  with  a  certain  power  of  tuQ  1  as 
shown  by  Minster  and  Vassiliou  [1979]  and  Minster  [1980]  (The  exponent  depends  on  the 
frequency  dependence  of  Q).  It  is  easy  to  verify  that  these  simple  scaling  laws  are  not 
obeyed  by  the  pulses  shown  here,  which  raises  the  question  of  how  the  high-strain 
laboratory  data  should  be  interpreted.  It  is  possible,  however,  that  scaling  laws  might  be 


11 


obeyed  reasonably  well  by  pulses  with  a  different  shape,  more  similar  to  those  generated  in 
the  laboratory  experiments. 

Discussion 

An  issue  of  some  importance  is  that  although  the  algorithms  satisfy  the  basic  laws  of 
continuum  mechanics,  it  is  not  guaranteed  that  they  simulate  the  rheology  of  any  real 
material.  The  application  of  amplitude-dependent  properties  in  a  completely  local  mode  — 
both  in  space  and  in  time —  remedies  objections  raised  to  the  treatment  of  the  same  problem 
by  Minster  and  Day  [1986].  However,  the  drop  in  displacement  pulse  spectral  amplitudes 
at  DC  is  somewhat  counter-intuitive,  intriguing,  and  peculiar  to  the  nonlinear  regime.  This 
behavior  is  quite  different  from  that  of  the  classical  linear  “ Q  operator”,  which  preserves  the 
area  under  a  unipolar  displacement  pulse  as  it  propagates  and  attenuates.  The  nonlinear 
calculation  does  not  preserve  this  area,  nor  does  it  preserve  unipolarity.  This  caused  us  to 
become  concerned  about  the  application  of  the  boundary  condition  at  one  end  of  the  one¬ 
dimensional  medium.  A  finite-duration  force  applied  to  a  boundary  would  immediately 
produces  a  displacement  pulse  with  a  tail  which,  strictly  speaking,  is  of  infinite  duration 
(i.e.,  decays  asymptotically  to  zero).  In  other  words,  we  introduce  a  finite-duration  pulse 
in  a  medium  that  is  incapable  of  supporting  one,  which  is  true  even  in  a  linearly  anelastic 
medium.  Thus,  our  finite-duration  displacement  pulse  implies  a  forcing  function  of  infinite 
duration,  and  rather  artificially  shaped  to  suppress  the  late-time  displacement  at  the  source. 
While  the  spectral  ratios  for  a  linear  Q  model  are  insensitive  to  the  nature  of  the  forcing,  in  a 
nonlinear  model  they  will  depend  on  its  time-history  and  amplitude.  There  exists  therefore 
a  possibility  that,  in  the  nonlinear  case,  the  boundary  condition  interacts  with  the  outgoing 
pulse  in  a  way  that  could  affect  the  spectral  analysis.  This  is  intimately  related  to  the  issue 
that  the  propagation  is  pulse-shape  dependent,  so  that  use  of  an  artificial  input  pulse  could 
cause  strange  results. 

As  a  check,  we  constructed  a  composite  one-dimensional  medium  by  adjoining  a  linear 
elastic  section  to  the  nonlinear  one.  A  displacement  pulse  of  the  form  (9)  was  then  created 
in  the  elastic  section,  and  the  total  linear  momentum  in  the  system  calculated  at  a  time  when 
the  trailing  edge  of  this  pulse  was  clearly  detached  from  the  boundary,  but  its  leading  edge 
had  not  yet  reached  the  nonlinear  section.  At  this  time,  no  more  work  was  performed  by 
the  boundary  condition,  and  the  end  node  could  be  left  free.  At  the  interface  between  the 
elastic  and  nonlinear  media,  this  incident  pulse  generates  both  a  reflected  and  a  transmitted 
pulse,  and  the  transmitted  pulse  was  attenuated  as  described  previously,  in  a  way  that 
preserved  the  total  momentum  in  the  system.  This  is  an  indication  that  the  simulations  are 
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dynamically  correct,  but  it  says  little  about  the  thermodynamics  of  the  system,  other  than  it 
is  clearly  lossy. 

At  first  glance,  a  somewhat  surprising  aspect  of  field  observations  is  that  linear 
attenuation  models  should  perform  reasonably  well  in  explaining  data  which  we  feel  certain 
to  be  contaminated  by  nonlinear  effects.  This  is  true,  for  example,  of  the  SALMON  data 
modeled  by  McCartor  and  Wortman  [1985,  1989],  of  the  COWBOY  data  used  by  Minster 
and  Day  [1986],  and  of  the  laboratory  data  for  pressed  salt  collected  by  Larson  [1982],  We 
feel  that  a  plausible  reason  is  that  the  main  effect  to  be  explained  is  the  combination  of  pulse 
broadening  and  amplitude  decay  with  increasing  distance:  linear  Q  models  offer  enough 
flexibility  to  accommodate  a  fairly  wide  range  of  observations,  provided  that  one  is  willing 
to  chose  rather  extreme  models  (e.g.  Q1  >0.1  in  the  nonlinear  zone).  On  the  other  hand, 
a  straightforward  spectral  analysis  — a  linear  technique—  of  SALMON  near-field  data  led 
McCartor  and  Wortman  [1989]  to  the  conclusion  that  effective  Q  is  anomalously  low  at  low 
frequencies,  in  qualitative  agreement  with  our  results  on  low-frequency  spectral  ratios. 
However,  since  such  Q  values  are  much  smaller  than  the  typical  laboratory  values 
measured  on  the  same  rocks,  the  physical  meaning  of  these  models  is  questionable.  In 
fact,  in  cases  of  extreme  attenuation,  it  is  debatable  whether  the  concept  of  Q  is  useful  in  the 
first  place.  In  the  case  of  data  collected  within  the  “elastic  radius” — such  as  the  SALMON 
data —  the  physical  explanation  offered  by  Sammis  [1989]  which  ascribes  the  pulse 
broadening  to  damage  mechanics,  as  developed  by  Ashby  and  Sammis  [1989],  seems  a 
good  candidate.  In  particular,  it  incorporates  scale  effects  explaining  why  comparable 
pulse-broadening  is  not  seen  in  small-scale  explosions  in  dry  rock  [e.g.  Nagy  and 
Florence ,  1986,  1987,  1988],  On  the  other  hand,  whether  this  explanation,  primarily 
intended  for  granite,  is  also  appropriate  for  salt  is  open  to  question. 

As  we  have  seen  (see  Figure  7),  the  puzzling  fact  that  Larson  [1982]  observed  apparent 
linear  superposability  of  pulses  in  salt  with  peak  strains  in  excess  of  10'4  is  not  so 
surprising  after  all.  This  superposability  test  is  simply  too  weak  a  test  of  linearity.  As 
argued  by  Savage  and  Hasegawa  [1967],  a  more  general  test  of  linear  superposition,  and 
thus  of  linearity,  would  be  achieved  by  comparing  the  attenuation  of  both  sinusoidal  and 
non-sinusoidal  disturbances  in  the  spectral  domain.  (They  used  this  technique  to  reject  the 
variable-friction  mechanism  originally  proposed  by  Knopojfand  MacDonald  [1958]  as  a 
candidate  for  attenuation  of  low  amplitude  seismic  waves.)  This  technique  might  be  adapted 
to  our  present  problem  by  comparing  spectral  decays  from  free  vibrations  of  a  test  bar  in 
the  laboratory  [e.g.  Bulau  et  al.,  1985],  with  those  estimated  from  travelling  pulses  as 
measured  in  Larson’s  [  1982]  experiments.  In  point  of  fact,  the  forced  resonance  curves 
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measured  by  Tittman  and  his  co-workers  in  essentially  the  same  material,  are  increasingly 
distorted  as  strains  increase  beyond  I  O'6,  a  clear  indication  of  nonlinearity. 

Implications  for  explosion  sources 

In  order  to  achieve  any  sort  of  predictive  capability  concerning  coupling  and  coupling 
variability  of  underground  explosions,  one  must  be  capable  of  constructing  numerical 
models  that  simulate  the  expected  radiation  field  under  various  hypotheses.  Considerable 
research  and  development  has  been  targeted  at  this  problem  over  the  past  30  years,  and 
models  have  reached  a  remarkable  level  of  sophistication  [see  for  example,  Cherry  et  al., 
1975;  Bache,  1982;  Cherry  and  Rimer,  1982;  Day  et  al.,  1983].  Most  of  this  effort  has 
been  aimed  at  modeling  accurately  the  regime  prevalent  immediately  around  the  cavity, 
where  extreme  nonlinear  damage  takes  place  because  the  loading  path  crosses  the  yield 
surface.  As  noted  earlier,  some  controversial  aspects  remain,  particularly  concerning  the 
effect  of  dilatancy  [e.g.  Nikolay evskiy  et  al.,  1978;  Cherry  and  Rimer,  1982;  Bache, 
1982;  Scholz,  1982],  The  potential  implications  of  the  simulations  described  in  the 
previous  section  for  the  characterization  of  explosion  sources  are  several.  We  shall 
summarize  them  in  this  section. 

If  we  assume  that  Y'l,3  scaling  with  yield  holds  [e.g.  Trulio,  1978;  Larson,  1982], 
nonlinear  rheology  may  well  prevail  to  scaled  slant  ranges  as  large  as  lft-IO3  mkg'113  [e.g. 
Minster  and  Day,  1986].  In  contrast,  numerical  models  of  explosive  sources  usually 
assume  the  rheology  to  be  linear  beyond  the  so-called  “elastic”  radius,  where  the  stress- 
strain  path  followed  by  material  particles  ceases  to  intersect  the  assumed  yield  curve;  this 
may  happen  at  strains  as  large  as  10~4  [e.g.  Cherry  et  al.  1975;  Cherry  and  Rimer,  1982]. 
Therefore,  there  has  to  exist  a  domain  surrounding  the  explosion,  where  the  peak  strains 
are  not  so  large  as  to  cause  failure  of  the  medium,  but  large  enough  to  trigger  nonlinear 
attenuation  mechanisms.  This  conclusion  was  reached  by  Minster  and  Day  [1986]  and  is 
not  invalidated  by  our  present  results. 

In  cavity-decoupling  scenarios,  complete  decoupling  is  defined  by  the  condition  that  the 
cavity  wall  is  not  subjected  to  stresses  high  enough  to  cause  it  to  yield  [e.g.  Larson, 
1985].  However,  this  definition  does  not  preclude  loads  high  enough  to  cause  nonlinear 
attenuation.  In  fact,  if  complete  decoupling  were  defined  in  terms  of  truly  linear  behavior 
of  the  cavity  wall,  then  it  w.mIJ  require  such  large  cavities  as  to  be  impractical  [e.g. 
Archambeau,  1985|.  These  issues  are  discussed  in  detail  in  the  proceedings  of  the  June, 
1985  DOE  workshop  held  in  Fajaro  Dunes,  California  [Larson,  1985].  An  important 
issue,  from  the  point  of  view  of  source  theory  is  whether  this  has  a  significant  effect  on 
observable  teleseismic  waves,  and  particularly  on  coupling  in  the  1  Hz  band,  since  this  is 
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what  controls  magnitude  estimates.  Some  of  the  results  presented  here  indicate  that  this 
might  indeed  be  the  case,  but  because  the  solution  to  one  nonlinear  problem  cannot  be 
safely  applied  to  another  problem,  there  remains  an  uncomfortable  margin  for  uncertainty. 

We  have  therefore  clear  reason  to  think  that  we  observe  the  seismic  source  through  a 
screen  consisting  of  a  nonlinear  zone  in  the  immediate  vicinity  of  the  source.  If  Q~l 
increases  very  strongly  with  increasing  amplitude,  (i.e.  the  medium  is  strongly  nonlinear), 
then  we  should  expect  amplitudes  to  decrease  very  rapidly  with  range  close  to  the  source, 
and  drop  to  linear  levels  almost  immediately.  On  the  other  hand,  in  the  case  of  mild 
nonlinearity,  high  amplitudes  will  persist  to  larger  ranges,  being  controlled  predominantly 
by  geometrical  spreading,  and  the  size  of  the  nonlinear  zone  could  be  much  larger.  In  the 
former  case,  intuition  would  indicate  that  the  dimension  of  the  nonlinear  attenuating  zone  is 
much  smaller  that  teleseismic  wavelengths,  so  that  teleseismic  waves  might  not  be  affected 
very  much;  in  the  latter  case,  the  converse  might  be  true.  Unfortunately,  the  disquieting 
result  shown  on  Figure  1 1  is  that  nonlinear  rheology  can  affect  the  spectral  content  of  a 
source  pulse  at  all  wavelengths.  This  means  that  even  a  small  nonlinear  zone — say,  a 
fraction  of  a  wavelength  thick —  has  the  potential  to  affect  significantly  the  effective  source 
spectrum  at  teleseismic  frequencies  of  interest.  In  particular,  we  cannot  rule  out  the 
possibility  that  body  wave  magnitudes  (and  therefore  yield  estimates)  might  be  affected  for 
all  explosions,  large  and  small.  To  complicate  matters,  if  the  relevant  physical  mechanism 
is  indeed  due  to  friction  between  crack  surfaces  in  the  rock,  then  one  should  have  to 
account  for  differences  in  attenuation  above  and  below  the  explosion,  due  to  the  differences 
in  confining  pressure;  this  suggests  that  if  nonlinear  attenuation  is  important  at  all, 
spherically  symmetric  models  will  not  suffice  to  describe  its  effects  on  the  far  field 
radiation. 

On  the  other  hand,  we  view  such  intuitive  arguments  with  considerable  caution  when 
nonlinearities  play  a  role.  In  particular,  we  must  emphasize  that  quantitative  results  obtain 
from  one-dimensional  simulations  cannot  be  used  to  infer  quantitative  conclusions  for  even 
the  simplest  (spherically  symmetric)  three-dimensional  case.  This  is  because  geometrical 
spreading  cannot  be  accounted  for  after  the  fact,  but  must  be  introduced  in  the  heart  of  the 
calculation,  where  it  does  not  commute  with  attenuation.  Furthermore,  we  have  noted  that 
the  spectral  ratio  results  shown  in  Figure  1 1  are  pulse-shape  dependent.  This  means  that, 
except  in  a  qualitative  sense,  firm  conclusions  concerning  the  impact  of  nonlinear 
attenuation  on  the  effective  source  will  require  true  three-dimensional  simulations  using  a 
realistic  source  function. 

One  constraint  on  any  nonlinear  rheological  model,  as  discussed  by  Minster  and  Day 
[  1986]  and  also  by  Trulio  [1985],  is  that  it  must  not  violate  the  empirically  verified  yield 
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scaling  laws.  Only  rate-independent  rheologies  will  not  destroy  scaling,  and  even  a  linear 
Q  model,  sensu  stricto ,  does  introduce  a  rate-dependence  through  the  physical  dispersion 
required  to  satisfy  causality.  This  means  that  large  explosions  which  generate  a  wavefield 
considerably  enriched  in  lew  frequencies  relative  to  small  ones,  “sec”  a  slightly  different 
medium  after  the  Y113  scaling  is  applied.  In  practice,  the  main  reason  why  a  linear  Q 
model,  even  with  extremely  low  Q  values,  does  not  invalidate  the  scaling  laws  is  that  the 
physical  dispersion  associated  with  attenuation  is  not  severe.  Since  our  calculations  show 
that  the  characteristics  of  nonlinear  pulse  propagation  are  pulse-shape  sensitive,  a  special 
set  of  three-dimensional  simulations  will  be  required  to  verify  the  scaling  is  preserved. 

On  the  other  hand,  Murphy  [1989,  personal  communication]  points  out  that  observed 
teleseismic  moment  estimates  for  nuclear  explosions  are  in  quite  good  agreement  with  those 
expected  on  the  basis  of  cavity  dimensions.  This  should  not  be  the  case  if  nonlinear  effects 
are  significant.  This  comment  has  considerable  merit,  but  is  not  necessarily  in  direct 
conflict  with  the  results  presented  above.  The  one-dimensional  simulations  do  not  deal  at 
all  with  near-field  effects,  and  only  show  that  the  final  (static)  displacement  field — 
including  the  final  cavity  size —  is  reached  according  to  a  history  (time  function)  that 
depends  on  range  in  a  more  complicated  way  than  we  typically  assume,  and  is  affected  by 
the  local  nonlinear  rheology.  Strictly  speaking,  the  moment  is  not  affected.  On  the  other 
side  of  the  issue,  Evernden  et  al.  [1986]  argue  that  the  equivalent  elastic  radii  inferred  from 
teleseismic  corner  frequency  measurements  for  many  explosions  in  different  media  are 
larger  than  those  predicted  from  finite  difference  models  adjusted  to  fit  near-field 
observations  le.g.  Rimer  and  Cherry,  1 982].  Further,  these  authors  used  a  very  simple 
source  representation  due  to  Sharpe  [1942],  in  which  a  step  function  pressure  of  amplitude 
Po  is  applied  at  an  elastic  radius  Rot  where  Po  should  be  comparable  to  the  effective 
compressional  strength  of  the  material.  They  conclude  that  values  of  Po  required  to  fit  the 
data  fall  well  below  the  rock  strengths  used  in  the  numerical  models,  and  the  values  of  Ro 
are  larger  than  the  numerical  “elastic”  radius.  This  may  suggest  either  unsuspected 
weakness  of  the  surrounding  material  near  explosions  detonated  at  NTS  and  other  U.S. 
sites,  or  unmodeled  nonlinear  behavior  extending  to  larger  distances  from  the  source. 

Conclusions 

Because  most  rocks  exhibit  nonlinear  amplitude-dependent  attenuation  at  moderate  to 
high  strains  and  relatively  low  confining  pressure,  the  transition  to  the  elastic  regime  in  the 
vicinity  of  an  underground  explosion  involves  a  nonlinear  zone  separating  the  classical 
“elastic  radius”  from  the  truly  linear  (anelastic)  far-field  regime.  Modeling  the  wavefield 
within  that  zone  is  difficult  because  we  cannot  use  the  usual  mathematical  techniques 
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which,  for  the  most  part,  depend  on  the  principle  of  superposition,  and  therefore  on  the 
assumption  of  linearity. 

The  simple  numerical  experiments  described  in  this  paper  lead  us  to  one  major 
conclusion,  namely  that  the  effects  of  nonlinear  attenuation  on  wave  propagation  are  not 
easily  predicted  by  simple  arguments,  even  after  simulation.  In  some  instances  of  interest, 
such  as  the  evaluation  of  Q~l  from  Lorentz  peaks,  nonlinearity  seems  to  have  only  a  modest 
effect.  Similarly,  some  departures  from  linear  pulse  superposition  are  not  very  severe.  On 
the  other  hand,  the  results  demonstrably  depend  on  the  initial  pulse  shape,  so  that  it  is 
necessary  to  assess  each  problem  through  a  separate  calculation  using  the  correct  geometry 
and  the  right  pulse  shape.  It  is  a  measure  of  our  lack  of  reliable  insight  into  this  class  of 
nonlinear  wave  propagation  problems  that  even  the  simplest  simulations  are  not  easily 
interpreted,  and  that  the  inferences  one  can  attempt  to  draw  concerning  more  realistic 
situauons  remain  somewhat  tentative. 
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Laboratory  Salt  Data 


Fig.  2.  Distorted  Lorentz  peaks  obtained  from  flexure  experiments  in  dome  salt  Data 
generated  by  Bulau  and  Tittman  [1983,  personal  communication].  Samples  were 
flexed  at  480  Hz  under  1.36  x  107  Pa  effective  pressure.  Each  resonance  curve 
has  been  rescaled  and  aligned  to  a  common  peak  for  this  plot,  in  order  to  illustrate 
the  distortion  at  large  strain  amplitude.  In  order  of  increasing  strain  amplitude, 
the  four  curves  correspond  to  a)  £  =  1.3  10~7,  Q  =  495,  b)  e  =  65  /O'7,  Q  =  403, 
c)  £  =  2.6  lO’6,  Q  =  319,  d)  £  =  6.3  10~6,  Q  =  242,  where  Q  is  estimated  from 
the  peak  width  at  half-power. 


Fig.  4.  Attenuation  as  measured  from  the  curves  shown  in  Fig.  3  using  the  relation  Q1 = 
((02-0)})/ (Or  at  the  half-power  level,  plotted  against  maximum  strain  amplitude. 
Apparent  Q'1  varies  linearly  with  emax-  The  slope  of  this  line  implies  a  /value  of 
about  150,  very  close  to  the  actual  value  of  140  used  to  generate  the  curves 
shown  in  Fig.  3. 
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Fig.  5.  Velocity  time  histories  at  four  propagation  distances,  for  a  nonlinear  simulation 
with  2cj  equal  to  1000  and  ?  x  equal  to  0.3.  A  plane  wave  is  incident  from 
one  direction,  as  shown.  The  prescribed  time  history  at  distance  0  is  shown  at  the 
far  left;  all  times  are  scaled  to  the  source  duration,  T,  and  all  distances  are  scaled 
to  the  product  of  wavespeed  and  source  duration,  cT. 
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Fig.  6.  Velocity  time  histories  at  three  propagation  distances,  for  a  nonlinear  simulation 
with  Qa  equal  to  1000  and  1  x  equal  to  0.3,  as  in  the  previous  figure.  In  this 
case,  however,  plane  waves  are  incident  from  both  directions.  Location  C  is  the 
midpoint  between  the  source  points.  Times  and  distances  are  scaled  as  in  the 
previous  figure. 
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Fig.  7.  (a)  The  velocity  time  histories  at  location  C  (midpoint).  Results  are  compared  for 

the  single  source  simulation,  multiplied  by  2  (dashed  curve),  and  the  interfering 
sources  simulation  (solid  curve),  (b)  Same  as  (a),  but  for  location  B. 


Fig.  8.  Displacement  time  histories  for  a  linear  (amplitude-independent  Q)  simulation, 
with  Q  equal  to  20.  The  source  time  history  is  the  derivative  of  that  shown  in 
Figure  5  (i.e.,  the  source  displacement  has  the  same  form  as  the  source  velocity 
shown  in  Figure  5). 
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Fig.  9.  Spectral  ratios  for  the  amplitude-independent  Q  simulation  shown  in  Figure  8. 
Each  curve  is  the  spectral  amplitude  of  the  corresponding  displacement  time 
history  from  Figure  8,  divided  by  the  spectral  amplitude  of  the  source 
displacement. 


31 


Time/Source  duration  (t/T) 


Fig.  10.  Displacement  time  histories  for  a  amplitude-dependent  Q  simulation,  with  (2a 
equal  to  1000  and  the  product  [?  x  e^]  equal  to  0.3.  The  source  prescription  is 
the  same  as  for  the  linear  simulation  depicted  in  Figure  8. 
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Fig.  1 1.  Spectral  ratios  for  the  amplitude-dependent  Q  simulations  shown  in  Figure  10 
Each  curve  is  the  spectral  amplitude  of  the  corresponding  displacement  time 
history  from  Figure  10,  divided  by  the  spectral  amplitude  of  the  source 
displacement. 
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1400  East  2nd  South 
Salt  Lake  City,  UT  84112 

Dr.  Stewart  W.  Smith 
Geophysics  AK-50 
University  of  Washington 
Seattle,  WA  98195 


Donald  L.  Springer 

Lawrence  Livermore  National  Laboratory 
L-205 

P.  O.  Box  808 
Livermore,  CA  94550 

D-.  George  Sutton 
Rondout  Associates 
P.O.  Box  224 
Stone  Ridge,  NY  12484 


Prof.  L.  Sykes 

Lamont-Doherty  Geological  Observatory 
of  Columbia  University 
Palisades,  NY  10964 


Prof.  Pradeep  Talwani 
Department  of  Geological  Sciences 
University  of  South  Carolina 
Columbia,  SC  29208 


Dr.  David  Taylor 
ENSCO,  Inc. 

445  Pineda  Court 
Melbourne,  FL  32940 


Dr.  Steven  R.  Taylor 

Lawrence  Livermore  National  Laboratory 

L-205 

P.  O.  Box  808 
Livermore,  CA  94550 

Professor  Ta-Liang  Teng 
Center  for  Earth  Sciences 
University  of  Southern  California 
University  Park 
Los  Angeles,  CA  90089-0741 

Dr.  R.B.  Tittmann 

Rockwell  International  Science  Center 

1049  Camino  Dos  Rios 

P.O.  Box  1085 

Thousand  Oaks,  CA  91360 

Dr.  Gregory  van  der  Vink 
IRIS,  Inc. 

1616  North  Fort  Myer  Drive 
Suite  1440 

Arlington,  VA  22209 

Professor  Daniel  Walker 
University  of  Hawaii 
Institute  of  Geophysics 
Honolulu,  HI  96822 


William  R.  Walter 
Seismological  Laboratory 
University  of  Nevada 
Reno,  NV  89557 


Dr.  Raymond  Willem  an 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 
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Dr.  Gregory  Wojcik 
Weidlinger  Associates 
4410  El  Camino  Real 
Suite  1 10 

Los  Altos,  CA  94022 

Dr.  Lorraine  Wolf 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 


Prof.  Francis  T.  Wu 
Department  of  Geological  Sciences 
State  University  of  New  York 
at  Binghamton 
Vestal,  NY  13901 

Dr.  Gregory  B.  Young 
ENSCO,  Inc. 

5400  Port  Royal  Road 
Springfield,  VA  22151-2388 


Dr.  Eileen  Vergino 

Lawrence  Livermore  National  Laboratory 
L-205 

P.  O.  Box  808 
Livermore,  CA  94550 

J.  J.  Zucca 

Lawrence  Livermore  National  Laboratory 
P.  O.  Box  808 
Livermore,  CA  94550 


GOVERNMENT 


Dr.  Ralph  Alewine  111 
DARPA/NMRO 
1400  Wilson  Boulevard 
Arlington,  VA  22209-2308 


Mr.  James  C.  Battis 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 


Dr.  Robert  Biandford 
A  FT  A  (71 T 

Center  for  Seismic  Studies 
1 300  North  17th  St.,  Suite  1450 
Arlington,  VA  22209-2308 

Eric  Chael 
Division  9241 
Sandia  Laboratory 
Albuquerque,  NM  87185 


Dr.  John  J.  Cipar 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 


Cecil  Davis 

Group  P-15,  Mail  Stop  D406 
P.O.  Box  1663 

I. -os  Alamos  National  Laboratory 
Los  Alamos,  NM  87544 


Mr.  Jeff  Duncan 
Office  of  Congressman  Markey 
2133  Rayburn  House  Bldg. 
Washington,  DC  20515 


Dr.  Jack  Evemden 
USGS  -  Earthquake  Studies 
345  Middlefield  Road 
Menlo  Park,  CA  94025 


Art  Frankei 
USGS 

922  National  Center 
Reston,  VA  22092 


Dr.  Dale  Glover 
DIA/DT-IB 

Washington,  DC  20301 
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Dr.  T.  Hanks 
USGS 

Nat’l  Earthquake  Research  Center 
345  Middlefield  Road 
Menlo  Park,  CA  94025 

Paul  Johnson 

ESS-4,  Mail  Stop  J979 

Los  Alamos  National  Laboratory 

Los  Alamos,  NM  87545 


Janet  Johnston 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 


Dr.  Katharine  Kadinsky-Cade 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 


Ms.  Ann  Kerr 
IGPP,  A-025 

Scripps  Institute  of  Oceanography 
University  of  California,  San  Diego 
La  Jolla,  CA  92093 

Dr.  Max  Koontz 
US  Dept  of  Energy/DP  5 
Forrestal  Building 
1000  Independence  Avenue 
Washington,  DC  20585 

Dr.  W.H.K.  Lee 

Office  of  Earthquakes,  Volcanoes, 
&  Engineering 
345  Middlefield  Road 
Menlo  Park,  CA  94025 

Dr.  William  Leith 
U.S.  Geological  Survey 
Mail  Stop  928 
Reston,  VA  22092 


Dr.  Richard  Lewis 

Director,  Earthquake  Engineering  &  Geophysics 
U.S.  Army  Corps  of  Engineers 
Box  631 

Vicksburg,  MS  39180 

James  F.  Lewkowicz 
GL/LWH 

Hanscom  AFU,  MA  01731-5000 


Mr.  Alfred  Lieberman 
ACDA/VI-OA'State  Department  Bldg 
Room  5726 
320  -21  st  Street,  NW 
Washington,  DC  20451 

Stephen  Mangino 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 


Dr.  Robert  Masse 
Box  25046,  Mail  Stop  967 
Denver  Federal  Center 
Denver,  CO  80225 


Art  McGarr 

U.S.  Geological  Survey,  MS-977 
345  Middlefield  Road 
Menlo  Park,  CA  94025 


Richard  Morrow 
A  CD  A/VI,  Room  5741 
320  21st  Street  N.W 
Washington,  DC  20451 


Dr.  Carl  Newton 

Los  Alamos  National  Laboratory 

P.O.  Box  1663 

Mail  Stop  C335,  Group  ESS-3 
Los  Alamos,  NM  87545 

Dr.  Kenneth  H.  Olsen 

Los  Alamos  Scientific  Laboratory 

P.  O.  Box  1663 

Mail  Stop  D-406 

Los  Alamos,  NM  87545 

Mr.  Chris  Paine 
Office  of  Senator  Kennedy 
SR  315 

United  States  Senate 
Washington,  DC  20510 

Colonel  Jerry  J.  Perrizo 
AFOSR/NP,  Building  410 
Bolling  AFB 

Washington,  DC  20332-6448 


Dr.  Frank  F.  Pilotte 

HQ  AFTAQTT 

Patrick  AFB,  FL  32925-6001 


Katie  Poley 
CIA-ACIS/TMC 
Room  4X16NHB 
Washington,  DC  20505 


Mr.  Jack  Rachlin 
U.S.  Geological  Survey 
Geology,  Rm  3  Cl 36 
Mail  Stop  928  National  Center 
Reston,  VA  22092 

Dr.  Robert  Reinke 
WL/NTESG 

Kirtland  AFB,  NM  871 17-6008 


Dr.  Byron  Ristvet 

HQ  DNA,  Nevada  Operations  Office 

Attn:  NVCG 

P.O.  Box  98539 

Las  Vegas,  NV  89193 

Dr.  George  Rothe 
HQ  AFTAC/TTR 
Patrick  AFB,  FL  32925-6001 


Dr.  Alan  S.  Ryall,  Jr. 
DARPA/NMRO 
1400  Wilson  Boulevard 
Arlington,  V A  22209-2308 


Dr.  Michael  Shore 
Defense  Nuclear  Agency/SPSS 
6801  Telegraph  Road 
Alexandria,  VA  22310 


Mr.  Charles  L.  Taylor 
GL/LWG 

Hanscom  AFB,  MA  01731-5000 


Dr.  Larry  Tumbull 
CIA-OSWR/NED 
Washington,  DC  20505 


Dr.  Thomas  Weaver 
Los  Alamos  National  Laboratory 
P.O.  Box  1663,  Mail  Stop  C335 
Los  Alamos,  NM  87545 


Glv'SULL 
Research  Library 

Hanscom  AFB  ,  MA  01731-5000  (2  copies) 


Defense  Intelligence  Agency 

Directorate  for  Scientific  &  Technical  Intelligence 

Attn:  DT1B 

Washington,  DC  20340-6158 


Secretary  of  the  Air  Force 
(SAFRD) 

Washington,  DC  20330 

AFTAC/CA 

(STINFO) 

Patrick  AFB,  FL  32925-6001 

Office  of  the  Secretary  Defense 

DDR  &  E 

Washington,  DC  20330 

TACTEC 

Battelle  Memorial  Institute 

505  King  Avenue 

Columbus,  OH  43201  (Final  Report  Only) 

HQ  DN'A 

Attn:  Technical  Library 

Washington,  DC  20305 

DARPA/RMO/RETRIEVAL 

1400  Wilson  Boulevard 

Arlington,  VA  22209 

DARPA/RMO/Security  Office 

1400  Wilson  Boulevard 

Arlington,  VA  22209 

Geophysics  Laboratory 

Attn:  XO 

Hanscom  AFB,  MA  01731-5000 

Geophysics  Laboratory 
Attn:  LW 

Hanscom  AFB,  MA  01731-5000 


DARPA/PM 

1400  Wilson  Boulevard 

Arlington,  VA  22209 


Defense  Technical  Information  Center 
Cameron  Station 

Alexandria,  VA  22314  (5  copies) 
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CONTRACTORS  (FOREIGN) 


Dr.  Ramon  Cabre,  S.J. 
Observatorio  San  Calixto 
Casilla  5939 
La  Paz,  Bolivia 


Prof.  Hans-Peter  Harjes 
Institute  for  Geophysik 
Ruhr  University/Bochum 
P.O.  Box  102148 
4630  Bochum  1,  FRG 

Prof.  Ey stein  Husebye 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 


Prof.  Brian  L.N.  Kennett 
Research  School  of  Earth  Sciences 
Institute  of  Advanced  Studies 
G.P.O.  Box  4 

Canberra  2601,  AUSTRALIA 

Dr.  Bernard  Massinon 

Societe  Radiomana 

27  rue  Claude  Bernard 

75005  Paris,  FRANCE  (2  Copies) 


Dr.  Pierre  Mecheler 
Societe  Radiomana 
27  rue  Claude  Bernard 
75005  Paris,  FRANCE 


Dr.  Svein  Mykkeltveit 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY  (3  copies) 
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FOREIGN  (OTHERS) 


Dr.  Peter  Basham 

Earth  Physics  Branch 

Geological  Survey  of  Canada 

1  Observatory  Crescent 

OM-m'H,  Oninrin  CANADA  K1 A  0Y3 

Dr.  Eduard  Berg 
Institute  of  Geophysics 
University  of  Hawaii 
Honolulu,  HI  96822 


Dr.  Michel  Bouchon 
l.R.l.G.M.-B.P.  68 
38402  St.  Martin  D'Heres 
Cedex,  PRANCE 


tar.  1  Elmar  Bungum 
NTNF/NORSAR 
P.O.  Box  5 1 

N-2(M7  Kjeller,  NORWAY 


Dr.  Michel  Campillo 
Ohservatoirc  de  Grenoble 
l.R.l.G.M.-B.P.  53 
38041  Grenoble,  FRANCE 


Dr.  Kin  Yip  Chun 
Geophysics  Division 
Physics  Department 
University  of  Toronto 
Ontario,  CANADA  M5S  1A7 

Dr.  Alan  Douglas 
Ministry  of  Defense 
Blacknest,  Brimpton 

Rending  RG7-4RS,  UNITED  KINGDOM 


Dr.  Roger  Hansen 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 


Dr.  Manfred  Henger 

Federal  Institute  for  Geosciences  &  Nat'l  Res. 
Postfach  510153 
D-3000  Hanover  51.  FRG 


Dr.  Fekadu  Kebede 

Geophysical  Observatory,  Science  Faculty 

Addis  Ababa  University 

P.O.  Box  1176 

Addis  Ababa.  ETHIOPIA 

Dr.  Tormod  Kvaema 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 


Dr.  Peter  Marshall 
Procurement  Executive 
Ministry  of  Defense 
Blacknest,  Brimpton 

Reading  FG7-4RS,  UNITED  KINGDOM 

Prof.  Ari  Ben-Menahem 
Department  of  Applied  Mathematics 
Weizman  Institute  of  Science 
Rehovot,  ISRAEL  951729 


Dr.  Robert  North 

Geophysics  Division 

Geological  Survey  of  Canada 

1  Observatory  Crescent 

Ottawa,  Ontario,  CANADA  K1A0Y3 

Dr.  Frode  Ringdal 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 


Dr.  Jorg  Schlittenhardt 

Federal  Institute  for  Geosciences  &  Nat'l  Res. 
Postfach  510153 

D-3000  Hannover  51,  FEDERAL  REPUBLIC  OF 
GERMANY 


Universita  Degli  Studi  Di  Trieste 
Facolta  Di  Ingegneria 

Istituto  Di  Miniere  E.  Geofisica  Applicata,  Trieste, 
ITALY 


Ms.  Eva  Johannisson 

Senior  Research  Officer 

National  Defense  Research  Inst. 

P.O.  Box  27322 

S-102  54  Stockholm,  SWEDEN 
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